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By taking into account the possibility of all the intralayer as well as the interlayer current order¬ 
ings, we derive an eight-band model for interacting electrons in bilayer graphene. With the numerical 
solution to the model, we show that only the current orderings between the same sublattice sites can 
exist within the range of the physical interacting strength. This result confirms our previous model 
of spin-polarized-current phase for the ground-state of interacting electrons in bilayer graphene that 
resolves a number of experimental puzzles. 
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I. INTRODUCTION II. FORMALISM 


Bilayer graphene (BLG) has been attracted much 
attention because of its potential application to new 
electronic devices—The experimental observations on 
high quality suspended BLG samples^— show that the 
ground state of the electrons at the charge neutrality 
point (CNP) is insulating with puzzling properties in 
the presence of external magnetic field: (i) the insulating 
gap can be closed by a perpendicular electric field of ei¬ 
ther polarity, 7 - (ii) the gap grows greatly with increasing 
magnetic field B about 46 times larger than the Zeeman 
splitting at B > 0.1 Tesla^ (hi) the state is particle-hole 
asymmetric in the presence of the magnetic field r' (iv) 
there is a peak structure in the electric conductivity at 
small B and at the CNP,— and (v) there are quantum 
Hall states at all integers from 0 to ±4— On theoret¬ 
ical aspect, a number of models for the ground state of 
the electrons in BLG have been proposed, which include 
the ferroelectric-layer asymmetric stated 1,2 or quantum 
valley Hall stated a layer-polarized antiferromagnetic 
state ; 14 ' 15 a quantum anomalous Hall stated - — a quan¬ 
tum spin Hall state , 16 ’ 18 a superconducting state in co¬ 
existence with antiferromagnetisnifi^ an ordered-current 
statej22r—, and a spin-polarized current state (SPCS) — - 
Besides these, the insulating gap has also been explained 
as electron scattering on boundaries between domains 
with different stacking (AB and BA)— - The earlier stud¬ 
ies based on the renormalization group approach pre¬ 
dicted the Lifshitz transition in the interacting electron 
system of BLG.— Amongst these models, only the SPCS 
can qualitatively explain the above experimental obser¬ 
vations. 

By the SPCS model^ there are spin-dependent cur¬ 
rents between the lattice sites due to the interactions 
between electrons. In our previous work, only the intra¬ 
sublattice current orders were taken into account. One 
then raises a question: are there any possible inter¬ 
sublattice current orderings? In this work, we will an¬ 
swer this question by deriving and solving an eight-band 
model that counts in all the possible current orderings of 
interacting electrons in the BLG. 


The Hamiltonian of the electrons is 

H = — ^ 'tjjcl^Cjcr + 77 y '' VijSnjSrij ( 1 ) 

ija ij 

where c] 0 .(Cj a ) creates (annihilates) an electron of spin a 
at site i, Uj is the hopping energy between sites i and 
j, 5rii = rii — n is the number deviation of electrons at 
site i from the average occupation n, and v’s are the in¬ 
teractions between electrons. We use the tight-binding 
model considering only the intra-layer NN [between a 
(a ') and b (6')] and interlayer NN (between b and a 1 ) 
electron hoppings given by t = 3 eV and t\ = 0.273 
eV, respectively . 26 ’ 27 We will use the mean-field approxi¬ 
mation (MFA) and therefore adopt effective interactions 
that include screening due to the electronic charge fluc¬ 
tuations. As in our previous model^ the effective inter¬ 
actions are given as 

= V(T) = r(l + 4.69,.r + ,^) (2) 

where r = \r\ with r as a vector from site i to site 
j, and q s = 27re 2 yo is the screening constant with 
Xo = ti ln4/-7r(aeo) 2 (and eo = y/3t/2) the polarizability 
by the random-phase-approximation (RPA). 29 The form 
of v(r) is consistent with the RPA in the limit r —> oo. 
In our previous work^ we have shown that except the 
intra-sublattice current ordering there are no charge and 
spin orderings with absence of the magnetic field. We 
here confine ourselves to the case of no external magnetic 
field and at the CNP. We therefore ignore the on-site in¬ 
teraction and simplify the discussion only on the current 
orderings. 

First, we explain why the current orderings may appear 
in the system. By the MFA, the interaction term in Eq. 
CD is approximated as 

2 ^VijSnMj ~ X w b'( c ^ c icr)4a c ^- ( 3 ) 

ij ija 
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FIG. 1: (Color online) Bond currents between atoms a (black) 
and b (white) in top layer of BLG. The purple lines denote 
the current from atom b to atom a, the red and green lines 
are respectively the left-turn and right-turn currents from a 
to b. The red and green circles with arrow indicate the cir¬ 
culations of the current in the hexagon bonds. The diamond 
area enclosed by the dashed lines is the unit cell. 


In general, the average (cia-c^) is not always real, but 
may contain an imaginary part^ 

( C-iaCju) — Rijcr T Hija ■ (4) 

The imaginary part I{j a here corresponds to a current, 
which is self-consistently determined by the MFA. The 
treatment of the current between the sites of same sub¬ 
lattice has been detailed in Ref. [ 22 ]. 

In addition to the intra-sublattice currents, we here 
consider the inter-sublattice currents between the sites 
in the same layer as well as in different layers. First, we 
analyze the currents between the two sublattices a and 
6, for example, on the top layer shown in Fig. 1. The 
consideration applies to the currents between other sub¬ 
lattices. In Fig. 1, we draw all the nearest-neighbor (NN) 
bond currents between the a (black) and b (white) atoms. 
The purple lines denote the current from atom b to atom 
a, the red and green lines are respectively the left-turn 
and right-turn currents from a to b. The hexagons con¬ 
sisting of green and purple bonds contain the right-turn 
current loops, while the hexagons of red and purple bonds 
bear the left-turn current loops. The current in the pur¬ 
ple bond equals to the sum of the currents in the red 
and green bonds because of the current continuity. In 
this process, the current density at each atom does not 
vanish, which is different from the case of the currents 
between the same sublattice where the current density 
vanishes. Since currents are in loop forms, the total cur¬ 
rent is zero. 

From Fig. 1, we see that the translational invariance of 
the system is broken. One then needs a superlattice with 


unit cell containing more a and b atoms to describe the 
electrons. We find that there are three a atoms and three 
b atoms in the unit cell as enclosed by the dashed lines 
in Fig. 1 for a single layer. In Fig. 2, we depict the unit 
cell for the BLG consisting of three a atoms ai, a 2 and 
as and three b atoms b\. 6 2 and 63 on the top layer and 
three a' atoms , a 2 and 03 and three b' atoms b \, b 2 and 
63 on the bottom layer. From the top view, the b atoms 
are just on top of the a' atoms, while the a atoms are 
at the hexagon center of the bottom atoms. The lattice 
constant defined as the distance between the NN same 
atoms (for example a atoms) is a ~ 2.4 A and interlayer 
distance d ~ 3.34 A. We will use the unit a = 1 for the 
length. We take the basis of the unit cell as 

ei = (3/2, ^3/2) 

e 2 = (-3/2, V3/2) 

where the numbers in the parenthesis are the x and y 
coordinates. The area of the unit cell is S = |ei x e 2 | = 
3-\/3/2. The basis of Brillouin zone (BZ) in the reciprocal 
space are given as 

vi = e 2 x z/S = (1/3,1/a/3) (5) 

v 2 = z x ei/S = (—1/3,1/V3) (6) 

where z is a unit vector perpendicularly pointing out of 
the plane. The reciprocal lattice points 2ni7rFi + 2n 2 7rv 2 
with n\ and n 2 integers and the first BZ are shown in Fig. 
3. Note that the Dirac points at the corners of the origi¬ 
nal first BZ as shown in Fig. 3 are now the lattice points 
of the reciprocal space. Therefore, the zero energy of 
noninteracting electrons now appears at the origin since 
the energy is a periodic function of momentum in the 
reciprocal space of the superlattice. 

We will still use the term, sublattice, to distinguish 
the atoms in the original BLG lattice. So, there are four 
sublattices a, b , a', and b'. However, for the superlattice, 
there are 12 atoms in each unit cell. The term ‘subspace’ 
will be used for describing a space consisting of one kind 
of these 12 atoms; there are 12 subspaces in the system. 


A. 12-band model 

For description of the problem in momentum space, we 
first define the operator in real space 

^ = (^ 1 ,^ 2 ) ( 7 ) 

with 

V'jo-l = ( a ]l^ b ]l^ a U^ b ]2^ a j3a^j3a) ( 8 ) 

^ja2 = ( a 'jl ( 9 ) 

where a'^ a (fcj^) creates an electron of spin <7 at p,th 
atom a ( b ) in jth unit cell. Then, the operator in mo¬ 
mentum space is given by 

i’ka = y= ^2 ^ ex P( - *£' J) ( 10 ) 
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FIG. 2 : (Color online) Unit cell of the bilayer graphene con¬ 
taining three a atoms (black) ai, 02 and 03 and three 6 atoms 
(white with red edge) bi, 62 and 63 on the top layer and the 
corresponding a' (red) and b' (blue) atoms on the bottom 
layer. The arrows denote the interlayer currents between the 
a and b' atoms. 


where N is total number of the unit cells in the super¬ 
lattice, and the momentum k is confined within the first 
BZ enclosed by the dashed line in Fig. 3. For the sake 
of description, we name an /rth atom of sublattice l in 
unit cell of the superlattice as the atom. For example, 
a k 2 a creates an electron of momentum k and spin a at 
a ,2 atom. 

In the momentum-space representation, we obtain the 
expression for the non-interacting Hamiltonian as 


Ho =J2^L H (fylpk* 

ka 

with 



/ (j— 0 0 \ 

C = —t\ [ 0 er_ 0 

\ 0 0 cr_/ 


( 11 ) 


( 12 ) 


and 


H s =-t 


0 

1 

0 

e i(fci+fe 2 ) 


0 

gi(fci+fc 2 ) 


1 0 e -i( k i+ k 2) 

0 1 0 

1 0 1 

0 1 0 

1 0 e~ ikl 

0 e ik2 0 


Q g—i(fcl+fc 2 ) 

1 0 

0 e~ ik2 

ifcl Q 

0 1 

1 0 


where cr_ is the Pauli matrix, Ay and A ’2 are the momen¬ 
tum components along v\ and V 2 , respectively. By the 
MFA for the interactions, the element of the self-energy 


is obtained as 


-'ll ,<T 

-‘[IIS 


(k) 


1 

N 




where c' kl {cki^a) creates (annihilates) an electron of 

momentum k and spin a at subspace, and v l ^(q) is 
the lattice Fourier component of the interaction, 

v%(<l) = ex P (-*9 • H) (I 3 ) 

R 

with r’^ v the vector from the l' u atom to the l^ atom in the 
unit cell and the i?-summation running over the positions 
of the unit cells in the superlattice. The element of the 
effective Hamiltonian is Hjj^l^k) + (k). Since there 

are 12 atoms in the unit cell of the superlattice, we thus 
have a 12-band model for the electrons. 

For the carrier concentration close to the CNP, only 
the energy levels close to zero momentum need to be 
considered. We here consider the system at the CNP. 
In our previous work, the elements of the self-energy at 
k = 0 have been proved to be pure imaginary numbers 
that come from the current orderings for the system at 



FIG. 3 : (Color online) Lattice points in the reciprocal space. 
The solid blue points are for the original honeycomb lattice 
with the area enclosed by the black solid line as the first Bril- 
louin zone (BZ). The red-edge circles are the additional points 
for the superlattice with unit cell containing three a-atoms 
and three 6-atoms. The two vectors represent the basis of the 
reciprocal lattice. The diamond area with dashed-line edge 
is the first BZ for the superlattice. The inner hexagon is the 
equivalent high symmetry first BZ. 
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the CNP. We therefore have 

k' 

~ (i4) 

k' 

where the second line comes from the fact that the av¬ 
erage ^m{c\., l , vrj Ck'ino) is sizable only when k! is close to 
0. To understand the physics of such an element of the 
self-energy, we consider the current 

I\iv = I m ( c ]'l'va C jl^) 

= Im [<4l ' V a C kl^) eiHj ~ j,) } 

k 

k 

where the same approximation as in Eq. d again has 
been used in the last line for not too large | j — j'\. But 
the equality in the last line is valid for the l and l[, 
atoms belonging to the same unit cell. Therefore, the 
self-energy can be expressed as 

£“>(0) = £“> « -iv« (0)J«>. (15) 

Note that Ijjf = 0. The remaining 6x11 current el¬ 
ements are not independent quantities, but satisfy the 
current continuity law and the BLG lattice symmetry. 
We analyze the current orders below. 

(1) Currents between the atoms in the same sublattice. 
For this case, we have obtained the result in our previous 
works , 22 i 23 

I^=I l>a sin(£-r*,) (16) 

where I\ = (47r/3,0) is a Dirac point, and I l,<7 is a con¬ 
stant satisfying I a a = —I b '’ a and I b ' a = —I a ' a . Since 
K is a lattice point in the reciprocal space, the current 
ijjf is a periodic function in the superlattice. It means 
that the current Ijjf is the same when changing the po¬ 
sition of the //tli atom to another unit cell but with vth 
atom fixed. Therefore, the formula given by Eq. csd 
represents not only the currents between the same type 
atoms within the unit cell but also between that in dif¬ 
ferent cells. Figure 4 shows the currents flow through 
an atom from/to the same type atoms on a layer. The 
current density in this case is zero at every atom. All the 
bond currents between the same type atoms constitute 
current loops in the lattice. Instead of considering I l ’ a 7 
we define the order parameter 

A l a = -V3v[ l 2 (0 (17) 

We then have 

= -£“4 = = iA l jV3. (18) 



FIG. 4: (Color online) Currents through the central atom 
(black) from/to the same type atoms on a layer. The cur¬ 
rents are incoming (outgoing) from (to) the white (red) atoms. 
There are no currents between the black and the blue sites. 


The relationship given by Eq. (fT51) can be obtained from 
Eqs. m and ©. 

(2) Currents between atoms in different sublattices. By 
viewing Fig. 1, we can obtain the formula for the bond 
currents between the sublattices a and b. The result can 
be extended to the case of currents between two atoms 
of any other two sublattices l and V except ll' = ab' and 
IV = ba !. For IV ^ ab' and ba ', we have 

I^ = I ll '-oo S (K v .f%), (19) 

where I 11 ,cr is a constant, and K v ' s are three Dirac points 
defined as 

A'i = (-271-/3,271-/^) 

K 2 = (4tt/3,0) 

K 3 = (27r/3, 27 t/v / 3). 

One can check that the currents given by Eq. m satisfy 
the continuity law. Again, the is a periodic function 
in the superlattice. In Fig. 5, some other bond currents 
between the atoms in different unit cells are also depicted, 
which give rise to the larger current loops. For IV = ab ', 
the bond currents between a and b' atoms are depicted 
in Fig. 2. We have 

7 a6V =/ a6V c os{ Q v . (20) 

where Qi = A), Q 2 = I\ 3 , and Q 3 = K 2 . For the 
currents between the sublattices b and a', since the two 
sublattices look like a same sublattice from top view, we 
obtain the formula as 

/^V=/6»Vsm ( ^-r^' ) . 


( 21 ) 
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FIG. 5: (Color online) The same as Fig. 1 but with some 
currents of longer bonds connected to the original atom de¬ 
picted. 


Define the order parameters 

A“' = -3t#(0)/“ V , lV?ba' (22) 

A bo ' = -x/3u 1 b f'(0)/ 1 6 “ V . (23) 

From Eqs. (USD, USD, (HD, and (USD, we can determine 
the relationship between the elements of the self-energy. 
For example, we have 


Af/3 


_ ^-\db,<r _ v-v(z6,f7 jc\ 

^11 — 2j 12 — — ^13 

^-\db,<T v- ~\db,(T jcy ^-\db,(7 

^21 ~ —^22 I L ~ ^23 

^-\db)(r jc\ _ ^-\db,(7 _ ^-\db,(7 

~^31 / Z — ^32 — ^33 


(24) 


and E b “ ,CT = —E“ b,CT . Here, we have used the fact that 
v“ b (0) = u“ b (0) since rfy appearing in the definition of 
u“ b (0) given by Eq. (fT3l) is a shift of a superlattice con¬ 
stant from T^ b or in addition then a rotation of angle 
±27 t/ 3; because the lattice is invariant under rotations 
±27 t/ 3, these operations on r^ b do not change the result 
of the r-summation. The consideration applies to other 
interactions. We therefore consider only the interactions 
V j2 v (0) with indices for the atoms b±, b[ and the a' 
atom just beneath b± as show in Fig. 2. 

All the above results can be summarized in the self¬ 
energy matrix 


E CT 


St s. 


tb 


J tb 


(25) 


trix due to the interlayer currents. E t is given by 



( 0 

A a6 

A“ 

A ab 

-A“ 

-2A“ b \ 


-K b 

0 

-K b 

A* 

2A ab 

-A b 

i 

-A“ 

A ab 

0 

-2A ab 

A “ 

A“ b 

3 

-K b 

"A* 

2A“ b 

0 


A* 


A % 

—2A“ b 

-A“ 

A ab 

0 

A ab 


\2A ab 

A b 

-A“ 6 

-A b 

_ A ab 

o / 


with A“ = -\/3A“ and A b = \/3A b . The expression for 
Ef, is obtained by the replacements A“ —> —A*, A b —>■ 
—A“, and A“ b —> A^ b ' in the above matrix E t . For E tb , 
we have 


y\aa 

A ab ' 

/ /\aa / 

-2A“ b ' 

_2A“ a ' 

A“ b ' ' 

0 

A f 

A ba ' 

A bb ' 

-A ba ' 

-2A bb ' 

Add' 

-2 A f 

-2A“ a ' 

A“ b ' 

Add' 

Af 

-A ba ' 

A bb ' 

0 

-2A bb ' 

A ba ' 

A bb ' 

-2A““' 

A“ b ' 

Add' 

^<7 

A“ 6 ' 

Add' 

^<7 

-2A ab ' 

A b /' 

-2A bb ' 

-A ba ' 

A bb ' 

0 

1 

*2 b 

<1 


where A ba ' = V3A b a a '. 


B. 8-band model 

To get an effective low energy Hamiltonian, we first 
expand the noninteracting Hamiltonian matrix H s at k = 
0 to linear k, 


H, = h° 


K 


with 


h u = - 



hi = 


—it 


(fci + fc 2 )cr_ 
(fci + k 2 )cr- 


— (ki + k 2 )cr+ - (fci + k 2 )a + 

0 k±(j— — k 2 a+ 

— kl<7+ + k 2 cr- 0 


The matrix h° can be written as h° = —tM ® ay with ® 
as the Kronecker product and M is a 3x3 matrix with 
all the elements = 1. The eigenvalues of M are 0 with 
dually degeneracy and 3. From the corresponding eigen¬ 
functions of these eigenvalues, we define the matrix 

i ( eia 1 A 

T = — e~ la e~ la 1 0 <r 0 

V 3 \ 1 e ia 1 J 


with a = 27t/ 3. With this matrix, we do transformation 
of H s 


T^H S T = e 0 


H 0 F 
ft D 


where E t and E b are self-energy matrices for the top and 
bottom layers, respectively, and E tb is the self-energy ma- 


(26) 
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with eo = a/ 3^/2, and 


H 0 

F 

D 


PxC 1 + Py (?2 0 

0 -PxCTi + p y cr 2 

-3 (p x - ipy)e~ ia a\ 

3 (p x + ip y )cr i 

—2v / 3(Ji -p y <J2, 


(27) 

(28) 
(29) 


The correction FZ) 1 FT is a second order in the momen¬ 
tum, and will be neglected. 

The Hamiltonian H (, k ) given by Eq. (fl2l) and the self¬ 
energy T, a given by Eq. (E>D are transformed by I 2 <8> T 
with I - 2 as the 2x2 unit matrix operating in the space 
of top and bottom layers. The matrix C in H ( k ) is un¬ 
changed under the transformation. The self-energy ma¬ 
trices are transformed as 


where p x and p y are the momentum components in the 
x — y orthogonal system 


Px = (&i - h)/3 

P y = (ki + k 2 )/Vs. 

Clearly, H 0 describes the energy bands close to zero with 
the two diagonal matrices expressing the low-energy non¬ 
interacting electrons in the two valleys, while D gives rise 
to the bands with energies about ±2-\/3eo far from the 
zero. The low energy band of D is fully occupied while 
the high energy band is empty. There is no contribution 
to the physical process from these two bands. These two 
bands will be ignored. The off-diagonal matrices F and 
Ft are high order corrections. To see it, we take further 
transformation of the right-hand side of Eq. (1261) using 

h 0 

-D-^t o-o 



/ -m t A fa 2 0 \ 

TtX t T = A fo 2 m t 0 , (30) 

\ 0 0 0 / 

( —m b A a J b 'a 2 0 \ 

T^T, b T = A^ b ' a 2 m b 0 , (31) 

V 0 0 0 / 

/-A b /a- - im a 0 \ 

T^ tb T = —im- a A b /a. 0 , (32) 

\ 0 0 0 / 


with 


m t 


m a 


(K o \ 

[o A b a ) ’ mb 

^aa ^ab 

0 A bb ' )■ 


V0 A"/ ’ 


with I 4 as the 4x4 unit matrix and obtain 


T 


H 0 F 
Ft D 


T\ = 


Ho 


FT^Ft 0\ 

0 DJ ' 


The zeros in the third row and third column of matrices 
in Eqs. stem from the current continuity law 

that results in zero of the total current. 


By keeping only the energy bands close to zero, we get an effective eight-band Hamiltonian as 


H e H( p ) = 


(PxG 1 +PyCT 2 - m t 

A ab o 2 

-(t 1 + A b f)a + 
\ im t 


A a J>a 2 -{h+A -im a/ \ 

-PxVi +p y cr 2 + m t -im- a ~{t\ - A b a a )ct_ 

im i _ a p x a i +p y a 2 - m b A^ b 'a 2 

~{t\ - A b a a ')a + A“ b 'a 2 -p x cri + p y a 2 + m b J 


(33) 


where we have set the unit of energy as eo = 1. Note that 
by setting all inter-sublattice current orders as zero, the 
formula of H e H can be transformed to a matrix with two 
4x4 matrices in the diagonal. The two 4x4 matrices rep¬ 
resent respectively the four-band Hamiltonian in the two 
valleys. Therefore, equation (l33l) reduces to our previous 
result. From Eq. rtddl) . it is clear that the current or¬ 
ders A“ b and A“ b give rise to the inter-valley couplings. 
This is because that the momentum of an electron in 
these current processes changes as its moving direction 
changes. We can neglect current order A^° since it is a 
correction much smaller than the hopping t\ between b 
and a'. 

With the effective Hamiltonian, we then find out the 


eigen-energy E\ (j>) and the wavefunction cj)^ (p) by solv¬ 
ing the equation 

H $ f (P)<t>v\{P) = £aG»)/va ip) (34) 

with pi, A = 1,2, • • • , 8. For briefty, we here suppress the 
spin index on the wave function and energy. With the 
result of the wavefunctions, the order parameters can be 
determined self-consistently. The formula for the order- 
parameters in terms of the wavefunctions are presented 
in Appendix A. 

For the current to be spin polarized, we have A l J = 
aA u . To satisfy this condition, we need to set A“ b = 
—A ab and A aa = A bb . Then under the transforms 
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a —> —a, a —> a', b —> and p y —>• —p y , the effective 

Hamiltonian is unchanged. 


III. NUMERICAL RESULTS AND DISCUSSION 

To numerically solve the model, we first need to deter¬ 
mine the seven interactions i>“ b ( 0 ) in the above equations 
for the order parameters. By considering the symmetry 
of the lattice, only three different values for the interac¬ 
tions are in question: u^( 0 ) = u“( 0 ), t;“ b ( 0 ) = 
and Vii (0) = v^\ (0) = v b \ (0). In Appendix B, we 
present the formula for these interactions in terms of 
sums over the original lattice points and show the re¬ 
lation between u“ 2 ( 0 ) and the interaction parameter in 
our previous workj^. In unit of eo = 1, the three interac¬ 
tions are determined as v™ (0) = 4.25, (0) = 4.77, and 

v \i (0) = 2.94. With (0) = 4.25 (so v s = 6.37 as given 
in previous work2£), we get A“ ~ 1 meV = Ao (that is 
the experimental data£) by setting all the inter-sublattice 
current orderings to zero. 

We have solved the model at various relative strengths 
v i\{^)/ v i 2 ( 0 ) an d v ii (®)/ v i 2 ( 0 ) and at zero temper¬ 
ature. To a wide range of the relative strengths 
u ii( 0 )/^i 2 ( 0 ) an d v^\ ( 0 )/v “2 ( 0 ), we find that the or¬ 
derings A“ b , A^ b ' , A“ a / and A bb ’ vanish. The inter¬ 
layer current ordering A“ b appears when the interac¬ 
tion v± b (0) is strong enough. Figure 6 shows the re¬ 
sult for A ab = A“ b of spin-up (er = +) electrons as a 
function of the strength v ab '/v aa = ( 0 )/vi% ( 0 ) with 

v ii(®)/ v i 2 (0) = 5- (Since here A“ b , A“ ,b / A“ a ' , and 
A bb vanish, the current orderings A“, A b , and A“ b can 
be spin polarized.) It is seen that A ab shows up as in¬ 
creasing the strength v ab /v aa over 1.33, while the orders 



ab '/ aa 
V /V 


FIG. 6: (Color online) Current order A“ as function of the 
relative interaction strength v ab /v aa with Wn (0)/v“2 (0) = 5. 
Inset: orders A“ and A b . 



FIG. 7: (Color online) Current orders A“ and A b as functions 
of the external potential it. 


A a = A“ and A b = A b + as shown in the inset in Fig. 6 
decrease sharply at this strength. The critical strength 
1.33 is larger than the physical strength 0.69 determined 
by the above parameters. At v ab /v aa < 1.33, only A a 
and A b are finite but all other orderings vanish. The re¬ 
sult for A a and A b actually coincides our previous one . 23 
From Fig. 6 , we conclude that (1) the intra-sublattice 
current orderings cannot coexist with the inter-sublattice 
current orderings, and ( 2 ) within the physical strength of 
the interactions, only the former can appear. 

We here discuss the effect of external electric field on 
the current orderings. Suppose an external electric field 
E is applied perpendicularly to the BLG plane. Then, 
there is a potential difference 2 u = Eed/e (with e a 
screening constant) between the two layers. For small 
it, the external held does not affect the appearance of 
A ab since the critical interaction v ab /v aa = 1.33 is very 
strong. On the other hand, the intra-sublattice current 
orderings delicately depend on the external held. Figure 
7 shows the intra-sublattice current orders A a and A b as 
functions of the external potential it. It is seen from Fig. 
7 that the two order parameters decrease with increasing 
the potential and eventually vanish. 

We have also studied the model dependence of the cur¬ 
rent orderings. The model is modified with including the 
so-called 73 hopping term between the a and b' atoms. 
Except the t and t\ hopping terms, the 73 term with 
73 ss 3ti/4 is the strongest one among the remaining 
hoppings 3 At low energy, this hopping gives rise to the 
elements = v 3 (p x + ip y ) ( K valley inter¬ 
layer coupling) and H^ f = = v 3 (-p x + ip y ) (■-K 

valley interlayer coupling) with v 3 = a/373/ 2 in the effec¬ 
tive Hamiltonian. By including this hopping, the critical 
interaction strength for the appearance of A ab is pushed 
to higher value v ab '/v aa ss 1.5. At v ab '/ v aa > 1 . 5 ; A ab/ 
is almost the same as that given in Fig. 6 , A a = 0 and 










FIG. 8: (Color online) Left: Functions f a {p) (main panel) 
and fb(p) (inset) at v ab /v aa = 1. Right: Function f a y(p ) at 
v ab '/v aa = 1.5. 

A b = 0. While at v ab '/v aa < 1.5, A ab ' = 0, A a and A 6 
are two finite constants with reductions (due to the addi¬ 
tional hopping on them) less than 2%. This result means 
the above conclusion obtained from Fig. 6 is unchanged. 

We end this section with checking the approximation 
in Eq. (fT4|) that the function Im(c£ ; , V(T cupa) is sizable at 
small k. To see it, we define the three possible nonvan¬ 
ishing functions 

fa{p) = 2\/3 J ^Im(4 1<T a fc2o -) 

fbip) = 2^ J ^lui(bl la b k2 a) 

fab'ip ) = 6 J {a{ la b' kla } 

where 9 is the angle of the momentum p. The behaviors 
of these functions are shown in Fig. 8 . The left panel 
presents the functions f a {p) (main panel) and fb(p) (in¬ 
set) at the interaction strength v ab /v aa = 1.. The right 
panel is f a b'{k) at v ab /v aa = 1.5. It is clear that these 
functions are nonvanishing at small p. This justifies our 
approximation in Eq. m- 

IV. SUMMARY 

By taking into account all possible current orderings, 
we have derived an eight-band model for interacting elec¬ 
trons in BLG. We have solved the model to a wide range 
of the interaction strength and found that among the 
inter-sublattice currents only the interlayer current be¬ 
tween a and b' atoms can appear when the interactions 
between electrons at these atoms are strong enough. The 
critical interaction strength for the current orderings be¬ 
tween the a and b' sublattices is about 2 times stronger 
than the physical interaction strength. We have thus 
confirmed that within the range of physical interaction 


strength only the intra-sublattice current orderings are 
possible. Physically, because the inter-sublattice cur¬ 
rent orderings break the translational invariance, it is 
not likely to appear at weakly interacting systems; much 
stronger strength of interactions is needed for construct¬ 
ing a more symmetry broken state. On the other hand, 
since the intra-sublattice current orderings maintain the 
system’s translational invariance, it can be realized at 
relatively weak interactions. The triangle lattices in the 
BLG are the favorite spaces for the current orderings. In 
such a space, the current orderings can survive in long 
range; since the currents flow along three dictions, the 
current density and the total current are zero satisfying 
the requirement of the equilibrium condition. 
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APPENDIX A: Order parameters 


Here, we express the formula for calculating the order 
parameters from the wavefunctions. Denote the eigen 
operator as < I>a(p). The original operator i/j k<7 can be 
expressed as 


il) kr7 = {h®T r )$(p)<S>(p) (Al) 


where T r is a 6x4 matrix obtained by deleting the last 
two columns in T, (j>{p) is an 8 x 8 matrix with elements 
</Va(p)j and $(p) is an 8 -component vector with com¬ 
ponents 4 >a(p). We here delete the last two columns in 
T because they are related in the transformation to the 
eigen-states of matrix D given by Eq. (1291) that has been 
dropped away in obtaining the eight-band model given 
by Eq. (1331) . Explicitly, the matrix is 
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With the transform given by Eq. GU), the order param¬ 
eters are then calculated as 


A“ = 

v?m 

2 N 

£[ 0 ia(?#ia(p) 

kX 

- </>3 A (p)^3a(p)]/a(p), 

A* = 

< 2 ( 0 ) 

2 N 

£[</ , 2a(p)^2a(p) 

- 4>lx(p)<l>ix(p)]h(p), 


kX 


A ab = 

vim 

N 

£ Im ^AW 2 A 

kX 

(p) + </>3a(p) < / > 4a(p) 


+ e~ 

m ^\x(p)m(p) + 

e “ < / > 3A(P) < / > 2A(p)]/A(p), 

K b ' = 

N £lm[4 A (p)<fc5 

x(p) + /4 A (p)08a(p) 


fcA 


+ e Za <f>t\(p)<f>8\(p) + e M 07 A (p)</> 6 A(p)]/A(p), 
A ““ = Vl1 ^ £ Iix #1a (P)^5A(P) + (t>3x(p)<t>7X (?) 

kX 

+ e ~ ioi 4>\\{p)fa\{p) + z ia <t>\ x ( p )^ 5 a ( p )]/ a ( p ), 
A bb' = V 1A°) ^Im[4 A (p)^>6A(p) + ^>4 a(p)^8a(p) 

kX 

+ e^^A^SAb) + e M ^ A (p)</> 6 A(p)]/A(p), 

A ab' = ( Q) £lm[</>:[ A (p)(/ 6 A (p) +^ a (p)^sa(p) 

kX 

+ e "* a ^lA (p)^ 8 a(p) + e M ^ A (p)</> 6 A(p)]/A(p), 


where / A (p) = /[E A (p)] = (<f>|(p)$ A (p)) is the Fermi 
distribution function. Here, the ^-summation can be re¬ 
placed with the p-integration according to 



with a cutoff p c = 1 for the largest momentum. We thus 
have obtained the self-consistent equations for determin¬ 
ing the order parameters. 


APPENDIX B: Interaction parameters 


We here derive the formula for these interactions so 
that they are given by sums over the original lattice 
points instead of over the superlattice ones. By so do¬ 
ing, we will clearly see how the present model reduces to 
the previous one. Denoting = p^„ + z with p^ the 
planar component and z = 0 or d in z direction respec¬ 
tively for intralayer or interlayer interactions, we expand 
v(R + r^ v ) as 


J 


v{\R + ^„\) = y^ v {q,z)exv[iq-(R + ^)\ 


= -j£ v (\Q n + z ) exp [i(Qn + q] ■ $v] exp[«<f'• R] 


V 


(Bl) 


q,n 


where V is the area of the graphene layer, R represents the position of an unit cell in the superlattice, ^-summation 
in the second line runs over the first BZ given by the smaller hexagon in Fig. 3, and n over all the points including 
the red-circles and the blue points in Fig. 3. Here, we have used the relation exp(*Qn • R) = 1. From the definition 
for v l p V {q) and Eq. (IBID , we get 

vf!l(q) = X] V (\R + ^uCl) exp(—ig • R) 

R 

= ^-£u(|<3„ + q\,z)exp[i(Q n + q) ■ p£ 

n 

with vq = a/ 3/2 as the area of the unit cell of the original honeycomb lattice. For we have 

<>) = ^-£ i v (Qm,z) + v(\Q m + K\,z)exp(iK ■ p^) +v(\Q m - K\,z)exp(-iK ■ ff v )\exp(iQ m ■ p^){B2) 

u m 

where m sums over the blue points in Fig. 3. Here, noting that exp(iQ m • p^) = 1 and cos(A' • p^„) = —1/2 (for 
p“ = ±1) for the atoms in the same sublattice, we have 


C(°) = 7^'52[ V (Qrn,z)-v{\Qm + K\,z)\, for = ±1. 


(B3) 
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As stated after Eq. (l24l) in Sec. II, we consider the vectors A^A s between the atoms in the small unit cell containing 
only four atoms ai, bi, a' (beneath bi) and b (. Then, for r being a position of the unit cell containing the four atoms, 
we have the expansion 

v{\r + ^ v \) = y^2v{\Q m + q\,z)exp[i(Q m + q)-^ v +iq-r\, for l £ V 

q,m 

1 - , 
v A) = y 2 ^v(\Qm + q\,z)exp(iq-f), for 1 = 1' 

q,m 

where the g-summation runs over the first BZ given by the larger hexagon in Fig. 3. By definition, we have the 
Fourier components 


rf'Av) = )eM-iq-r) 

r 

= ^-'52v(\Qm + q\,z)exp[i(Q m + q)-ftl], for l ± l' (B4) 

= ^2v(r)exp(-iq-r) = ^-^2v(\Qm+q\,z). (B5) 

r ^ m 


By applying Eq. (IB4I) in Eq. (1B2[) , we get for l ^ l' 

«£( 0 ) = livU°)+<AK) + vU-K)] 

= i^H|r + ^:i)[ 1 + 2c °s(^-r)]. (B6) 

r 

From Eqs. (1B3I) and (1B5I) . we have 

«i‘ 2 (0) = §[*&(0 )-v l lAK)] 

= lY^vm-coAK-rj] 

r 

= |y^w(r)sin 2 3 4 5 6 * 8 (A'-f) (B7) 


where in the last line we have used the fact cos (A' • r) = 
cos(2 K ■ r) since the difference between 2 K and —K is a 
lattice vector (blue) in the reciprocal space. Using the in¬ 
teraction parameter v s as defined in our previous work^ 
we have tq 2 (0) = 2v s /3. Applying this result in the ex¬ 
pressions for A“ and Aj( in Appendix A and recognizing 

EJa(p)[4 2 )a(p)^(2)aW - 4(4)a(p)^3(4)a(p)] as the 
valley difference of the electron distribution function in 
a (6) sublattice, we reach the same forms for A“ and A^. 
as in our previous work^ 
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